The shape and erosion of pebbles 
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The shapes of flat pebbles may be characterized in terms of the statistical distribution of curvatures 
measured along their contours. We illustrate this new method for clay pebbles eroded in a controlled 
laboratory apparatus, and also for naturally-occurring rip-up clasts formed and eroded in the Mont 
St. -Michel bay. We find that the curvature distribution allows finer discrimination than traditional 
measures of aspect ratios. Furthermore, it connects to the microscopic action of erosion processes 
that are typically faster at protruding regions of high curvature. We discuss in detail how the 
curvature may be reliable deduced from digital photographs. 

PACS numbers: 45.70.-n,83.80.Nb,91.60.-x,02.60. Jh 
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I. INTRODUCTION 

The roundness of pebbles on a beach has long been a 
source of wonder and astonishment for scientists in many 
fields pJQ Explanations for the pebble shapes were born 
from the simple pleasure of understanding nature but also 
from the hope that a pebble, or a collection of pebbles, 
might carry lithographically imprinted the signature of 
their erosion history. Reading that imprint would then, 
for instance, reveal if a pebble was eroded on a beach, a 
river or a glacier, or if it traveled a long distance down 
a stream. It even perhaps would reveal for how long 
the erosion forces have been at work on that object. Of 
obvious interest in Geology Q, a physical understand- 
ing of the formation of erosion shapes would also allow 
for a better control of many industrial processes lead- 
ing to rounded objects such as gem stone or clay bead 
grinding in tumblers or fruit and vegetable peeling in 
several mechanical devices. Diverse mathematical tools 
have been developed for geometrical shape analysis of 
crystalites, cell membranes, and other far from equilib- 
rium systems 0, IE IS 0; however, these do not seem 
applicable to pebbles. 

The evolution of a pebble shape under erosion can ar- 
guably be viewed as a succession of elementary cuts that 
act at the surface of the body to remove a given amount 
of material. This converts young, polyhedral-like shapes 
with a relatively small number of large sides and sharp 
vertices into more mature shapes with a high number of 
small sides and smooth vertices. The size and the shapes 
of each of these successive ablations, as well as the surface 
sites where the cutting happens, are determined both by 
the conditions under which erosion takes place and by 
the nature of the material being eroded. Exposure of a 
young, polyhedral-like shape to the rough tumbling of a 
steep stream slope will result in relatively large cuts of 
the angular sections, while exposure to the gentle erosion 
of wind or water is more likely to lead to small cuts al- 
most parallel to the existing flat sides. Also, the same 
sequence of external forces acting on two identical origi- 
nal shapes of different materials will result into distinct 



forms due to weight, hardness or anisotropy differences. 
In spite of the diversity of factors at play in shape mod- 
ification, the complete evolution of the pebble shape is 
fully determined by (i) the initial form described by some 
number of faces, edges and vertices and (ii) the position, 
size and orientation of the successive ablations. 



Given that the erosion process evolves by a succession 
of localized events on the pebble surface, it is surprising 
that the majority of the precedent attempts to charac- 
terize the pebble shapes were restricted to the determi- 
nation of global quantities such as the pebble mass or 
the lengths of its three main axes Clearly, in order 
to capture both the local nature of the erosion process 
and the statistical character of the successive elementary 
cuts, one needs to build a new detailed description of the 
pebble shapes based on quantities that are more micro- 
scopic and more closely connected to evolution processes. 
In Ref. |S| we proposed curvature as a key microscopic 
variable, since, intuitively, protruding regions with large 
curvature erode faster than flatter regions of small cur- 
vature. We then proposed the distribution of curvature 
around a flat, two-dimensional, pebble as a new statisti- 
cal tool for shape description. And finally we illustrated 
and tested these ideas by measuring and modeling the 
erosion of clay pebbles in a controlled laboratory appa- 
ratus. 



In this paper we elaborate on our initial Letter 
and we apply our methods to naturally-occurring rip- 
up clasts found in the tidal flats of the Mont St. -Michel 
bay. Section II begins with a survey of shape quantifi- 
cation for two-dimensional objects, in general, and reca- 
pitulates our new curvature-based method. Section III 
provides further details of the laboratory experiments on 
clay pebbles. Section IV presents a new field study of the 
Mont St. -Michel rip-up clasts. And finally, following the 
conclusion, two methods are presented in the Appendix 
for reliably extracting the local curvature from digital 
photographs. 
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II. 2D SHAPE QUANTIFICATION 

The issue of rock shape is of long-stan ding interest in 
the field of sedimentology H El El El El El El El 

. Two basic methods have become sufficiently well es- 
tablished as to be discussed in introductory textbooks Q . 
The simplest is a visual chart for comparing a given rock 
against a standard sequence of rocks that vary in their 
sphericity and angularity. A rock has high "sphericity" 
if its three dimensions are nearly equal. And it is "very 
angular", independent of its sphericity, if the surface has 
cusps or sharp ridges; the opposite of "very angular" is 
"well rounded". While useful for exposition, such verbal 
distinctions are subjective and irreproducible. The sec- 
ond method is to form dimensionless shape indices based 
on the lengths of three orthogonal axes. From the ratios, 
and the ratios of differences, of the long to intermediate 
to short axes, one can readily distinguish rods from discs 
from spheres. And a given rock may be represented by a 
point on a triangular diagram according to the values of 
three such indices, with rod / disc / sphere attained at 
the corners. This practice is nearly half a century old [l8| . 
Nevertheless, there is still much debate about which of 
the infinite number of possible shape indices are most use- 
ful |19l l2Cl l2lL l22l |23| . In any case, such indices cannot 
capture fine distinctions in shape, let alone the verbal 
distinctions of angular vs rounded. Furthermore, they 
provide no natural connection to the underlying physical 
process by which the rock was formed. 

More recent methods of shape analysis employ 
Fourier [2 4 l25l 12(1 |2jj , or even wavelet |28| , transforms 
of the contour. This applies naturally to flat pebbles 
or grains, but also to flat images of three-dimensional 
objects. The advantage of Fourier analysis over shape 
indices is that, with enough terms in the series, the exact 
pebble contour can be reproduced. For simple shapes, 
the contour may be described in polar coordinates by ra- 
dius (eg distance from center of mass) vs angle, r(9), and 
the corresponding transform. However, this representa- 
tion is not single-valued for complex shapes with pits 
or overhangs. Generally, the contour may be described 
by Cartesian coordinates vs arclength, {x(s),y(s)}, and 
the corresponding transforms. In any case, the relative 
amplitudes of different harmonics give an indication of 
shape in terms of roughness at different length scales. In 
a different area of science, Fourier representations have 
proven especially useful for analysis of fluctuations and 
instabilities of liquid interfaces, membranes, etc. 0,H@- 
In practice, for these systems, shape fluctuations are 
sampled during some time interval and then the aver- 
age Fourier amplitudes extracted by averaging over many 
different realizations of the shape. Also, because these 
phenomena are linear, each Fourier component grows or 
shrinks at some amplitude-independent rate and the evo- 
lution is fully determined by a dispersion relation. Un- 
fortunately these features do not hold for the erosion of 
pebbles. Because each pebble shape only provides one 
configuration, average quantities need to be built from 



a different prescription. Also, there is no a priori guar- 
anty that the variables are Gaussian distributed, and one 
needs a direct space method to better assess the impor- 
tance of non-linear phenomena. Non-linearity is, we sus- 
pect, intrinsically embedded in the erosion mechanisms 
of pebbles. If one considers for instance a shape repre- 
sented by a single harmonic in the r{6) representation, 
it is clear that the peaks will wear more rapidly than 
the valleys. Therefore the erosion rate cannot be a func- 
tion of the harmonic number only; it must either be a 
non-linear function of the amplitude itself or a function 
coupling many harmonics. 

Our aim is to provide an alternative measure of pebble 
shape that is well-defined, simple, and connects naturally 
to local properties involved in the evolution process. We 
restrict our attention to flat pebbles, where an obvious 
shape index is the aspect ratio of long to short axes. 
Since erosion processes generally act most strongly on 
the rough, pointed portions of a rock, we will focus on 
the local curvature of the pebble contour. Technically, 
curvature is a vector given by K = dT/ds, the deriva- 
tive of the unit tangent vector with respect to arclength 
along the contour |29j . More intuititively, the magnitude 
of the curvature is the reciprocal of the radius of a circle 
that mimics the local behavior of the contour. Here we 
shall adopt the sign convention K > where the contour 
is convex (as at the tip of a bump) and K < where 
the contour is concave (as where a chip or bite has been 
removed from an otherwise round pebble). In the Ap- 
pendix, we describe two means by which the curvature 
may be reliably measured at each point along the peb- 
ble contour. Note that the average curvature is simply 
related to the perimeter of the contour: 

P = 2n/{K), (1) 

which is obviously correct when the shape is a circle. 

To describe the shape of a pebble, a very natural quan- 
tity is the distribution of curvatures, p(K), defined such 
that p(K)dK is the probability that the curvature at 
some point along the contour lies between K and K + dK 
8]. In order to distinguish different distributions, as a 
practical matter, it is more reliable |30j to use the cumu- 
lative distribution of curvatures 

f(K) = f K p(K')dK' (2) 
Jo 

Literally, f(K) is the fraction of the perimeter with cur- 
vature less than K. Note that f(K) increases from 
to 1 as K : — > 00; the minimum curvature is where 
f(K) first rises above 0, the maximum curvature is where 
f(K) first reaches 1, and the median curvature is where 
f(K) = 1/2. Unlike for p(K), it is not necessary to bin 
the curvature data in order to deduce f(K). Instead, just 
sort the curvature data from smallest to largest and keep 
a running sum of the arclength segments, normalized by 
perimeter. Finally, so that the shapes of pebbles of dif- 
ferent sizes may be compared, it is useful to remove the 
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FIG. 1: (a) Radius, (b) normalized curvature and (c) frac- 
tion of the perimeter f(K) with curvature less than K, vs 
K divided by the average curvature (A'), for a superellipse, 
oval, ellipse, and circle. The curve types match those for the 
shapes as shown in the inset. Note that, except for the circle, 
all have the same aspect ratio. The differences in shape are 
reflected in differences in the forms of f(K). 



scale factor (K), which is related to the total perimeter 
as noted above in Eq. (JIJ. Altogether, we thus propose to 
quantify pebble shape by examining f(K) as a function 
oiK/(K). 

To help build intuition, examples of f(K) are given in 
Fig-fflfor a few simple shapes. The simplest of all is a cir- 
cle, where the curvature is the same at each point along 
the contour: K = (K) = 2tt/P = 1/R. Thus f(K) = 
0(1) for K < (>)(K). The curvature distribution is the 
derivative of this step function, giving p(K) — 5(K—(K) ) 
as required. The other three shapes shown are a superel- 
lipse, an oval, and an ellipse. For the latter, with long and 
short axes a and b respectively, one may compute f(K) = 
E{sm~ 1 [l/e{ l-(2VT^ E(e 2 )/(TTK)) 2 / 3 ) 1 / 2 },e 2 )/E(e 2 ) 
where e = y/l — b 2 /a 2 is the ellipticity, E(x) is the com- 
plete elliptical integral of the first kind and E(x,m) is 
the incomplete elliptic integral of the second kind 29]. 

While noticeably different, the superellipse, ellipse, 
and oval in Fig. ^ all have the same aspect ratio of 1.5. 
This emphasizes how a single number is insufficient to 
quantify shape. The shape differences do, however, show 
up nicely in the forms of f(K). The ellipse is closest to 



a circle, with a distribution of curvatures that is most 
narrowly distributed around the average and hence with 
an f(K) that is most like a step function. The superel- 
lipse is farthest from a cirlce, with four long nearly-flat 
sections and four high-curvature corners; its curvature 
distribution is broadest. The oval is intermediate. 



III. LABORATORY EXPERIMENT 

The examples given in Fig. ^ correspond to regular, 
highly symetric shapes of two dimensional convex "peb- 
bles" . In practice, natural or artificial erosion processes 
lead to curvature functions with an important statistical 
component. In this section we elaborate on the labora- 
tory experiments of Ref . Q , designed to study both the 
statistical nature of the curvature distribution and the in- 
fluence of the original shapes of the pebbles on the final 
output of a controlled erosion process. 

Laboratory pebbles were formed from "chamotte" clay, 
a kind of clay made from Kaolin and purchased from 
Graphigro, France. The water content of the purchased 
clay was 22% in a state that could easily be kneaded. The 
clay was kept tightly packed before use in order to avoid 
water evaporation. Clay pebbles were produced using 
aluminium molds made in our laboratory. They consist 
of a polygonal well of 0.5 cm depth. Once the mould 
was filled with clay, it was left at rest for 24 hours, so 
that 98.5% of the water was removed by evaporation. 
All the experimental results presented here concern one 
day old pebbles. We noticed that pebbles older than 2 
days were too fragile for our experimental configuration. 
The dimensions of the various pebbles studied here are as 
follows: 4 squares of side 5 cm, 5 rectangles of sides 4 x 
6 cm 2 , 5 regular pentagons of side 4.25 cm, one triangle 
with sides 7 cm, 7. 5 cm and 9 cm one irregular polygon 
with 7 sides, one lozenge of acute angles 45° and side 
5 cm and one circle of diameter 7 cm. 

The wearing method that we chose relies on placing a 
pebble in the rotating apparatus sketched in Fig.|2 The 
apparatus is a square basin, of dimensions 30 x 30 x 7 cm 3 . 
The basin bottom is a 1 cm thick aluminium plate and 
the walls are made of 0.04 cm thick aluminium sheets. 
This rotating plate is fixed to a rod held by the jaws 
of a laboratory mixer, Heidolph RZR1. The mixer itself 
is fixed to a tripod, so that its inclination angle can be 
varied. 

A typical trajectory of the pebble during the contin- 
uous rotation of the plate can be described as follows. 
First the pebble rotates with the basin until it reaches 
a high position. After the plate has rotated an angle 
between 7r/2 and ir, the pebble begins to slide due to 
gravity, until it hits one of the walls in the bottom part 
of the container, and then rolls down along that wall as 
the basin keeps its rotation. After a short stop at a con- 
tainer corner the pebble starts a new cycle again. We 
performed preliminary tests in order to determine both 
ideal basin orientation and ideal rotation frequency for 
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FIG. 2: The wearing apparatus used for the laboratory ex- 
periments. The rotating metal tray is 30 x 30 x 7 cm 3 . 



our experiments. As expected, above some maximum ro- 
tation frequency the pebble becomes immobilised in the 
basin: centrifugal forces maintain the pebble on a given 
position against the wall. Also, under some minimum 
inclination angle, no fall of the pebble is observed, while 
a high inclination doesn't allow the pebble to reach it's 
maximum altitude. Altogether, we found it suitable to 
operate at a basin angle of 45° and a rotation frequency 
of one cycle per second. Using the latter experimental 
conditions, we observed that the in-plane dimension of 
a pebble decreased by around a factor 2 after 30 min- 
utes. Thus, a significant wearing of a pebble could be 
observed after a few minutes rotation. In practice, each 
pebble was eroded under the described conditions during 
30 minutes, while a picture of the pebble was taken after 
each 5 minutes wearing. Hence, for each of the pebbles 
studied, we obtained about 7 pictures, corresponding to 
the initial pebble and to six following states of the wear- 
ing pebble. For some of the pebbles 8 or 9 pictures at 5 
minutes interval were taken. The images were then ana- 
lyzed following the method described in the appendix. 

An example for the shape evolution produced by this 
method is given in Fig. El with photographs shown every 
five minutes. The corresponding cumulative curvature 
distributions f(K) are given in Fig. 0] where the inset 
shows the extracted contours. Here the initial shape 
is square, with four long nearly-flat regions and four 
short high-curvature regions. Thus the initial f(K) rises 
steeply around K = and extends with relatively little 
weight out to K ^> (K). At first, the action of ero- 




FIG. 3: Shape evolution of a 5 x 5 cm square pebble eroded 
in the laboratory, by the method explained in the text. 
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FIG. 4: Cumulative curvature distribution, f(K), for the 
evolving pebble depicted in the inset (and pictured in Fig. 3). 
As the pebble becomes progressively rounder, the curvature 
distribution narrows and approaches a final shape. The time 
interval between successive contours is 5 minutes. For the 
inset, the axes are given in pixel units, equal to 0.132 mm. 
For contrast, a circle is shown by points. 



sion is most rapid at the high-curvature corners, with 
the flat regions in between relatively unaffected. Thus 
the high- AT tail of f(K) at first is suppressed, and weight 
builds up across (0.5 — 2)(K). Next the rounded cor- 
ners erode further and gradually extend across the flat 
sections. Thus weight in f(K) is gradually concentrated 
more and more toward (K). After about 15-20 minutes, 
when the flat sections are nearly gone, the form of f(K) 
fluctuates slightly but ceases to change in any system- 
atic manner. In otherwords, the shape of the pebble has 
reached a final limiting form. Further erosion will affect 
pebble size, but not pebble shape! 

To test the universality of the final shape, we repeat 
the same experiment both for other squares as well as for 




FIG. 5: The initial and final forms of different shapes eroded 
in our experiment. Final shapes for all cases are not only sim- 
ilar in shape and in size but they have also the same curvature 
distribution. 



a variety of other initial shapes such as rectangles, trian- 
gles, and circles. A number of different examples showing 
both the initial and final shapes are shown in Fig. El In all 
cases, the cumulative curvature distribution f(K) shows 
a systematic evolution at short times and slight fluctua- 
tions about some average shape at later times, just as in 
Fig- El The more angular or oblong the initial shape, the 
more erosion is needed to reach a stationary final shape. 
The average final f(K) is shown for the various initial 
shapes in Fig. |U Evidently, these all display the same 
quantitative form independent of the initial shape. Even 
f(K) for an initially-circular pebble broadens from a step 
function to the same form as all the others. 

The final f(K) for all initial shapes can thus be av- 
eraged together for a more accurate description of the 
stationary shape produced by the laboratory erosion ma- 
chine. The result is shown by the open circles in the same 
plot, Fig. El Differentiating, we obtain the actual curva- 
ture distribution, p(K), shown on the right axis, ft is 
fairly broad, with a full-width at half-maximum equal to 
about 1.6(K). The actual shape is not quite symmetrical, 
skewed toward higher curvatures. The closest simple ana- 
lytic form would be a Gaussian, exp[— (K— (K)) 2 /(2a 2 )}. 
The actual distribution is slightly skewed toward higher 
curvatures, but the best fit gives a standard deviation of 
a = Q.70(K), as shown in Fig. It is easy to imagine 



FIG. 6: Integrated curvature distribution, f{K), for the final 
shapes of laboratory pebbles with various initial shapes, as 
labeled. These are indistinguishable to within measurement 
accuracy; their average is shown by the open circles. The cor- 
responding average curvature distribution, p(K) is obtained 
by differentiation and is shown on the right axis along with a 
fit to a Gaussian shape. 



that the width of this distribution could be set by the 
strength of the erosion process. For example, if the angle 
of the rotating pan were lowered, then the erosion would 
be more gradual and more like polishing; in which case 
a rounder stationary shape may be attained with a nar- 
rower distribution of curvatures. The form of f(K), as 
well as its width, could also be affected. These types of 
questions can be addressed, both in laboratory and field 
studies, now that we have an incisive tool like f(K) for 
quantifying shape. 

To further study the erosion produced by our labo- 
ratory apparatus, we now consider how the perimeter 
of the pebble decreases with time, P(t). Since the ini- 
tial behavior depends on the specific initial shape, we 
focus on subsequent erosion once the (universal) station- 
ary shape is achieved. If the final stationary shape of 
the curvature distribution is reached at time to, then 
the quantity of interest is really P(t) / P(to) vs t — to- 
The results, averaged over all laboratory pebbles, are 
shown in Fig. \7\ Though the dynamic range is not 
great, the data are consistent with an exponential de- 
crease, P(t) = P(t )exp[-(t - to)/r}. The best fit to 
this form is shown by a solid curve; it gives a decay con- 
stant of r = 44 min. Exponential erosion is, in fact, 
observed in field and laboratory studies [l^. It is to be 
expected whenever the strength of the erosion is propor- 
tional to the pebble size, as in our lab experiments where 
the impulse upon collision is proportional to the pebble's 
weight. 



IV. FIELD STUDY 

As a first field test of our method of analysis, we col- 
lected mud pebbles in the Mont St. -Michel bay, France. 
The littoral environment located at the inner part of the 
Norman-Breton Gulf is characterized by a macro-tidal 
dynamics. This location exhibits the second largest tide 
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FIG. 7: Perimeter vs time, where to is when the stationary 
shape has been reached. 




FIG. 8: Typical shapes of Mont St. -Michel rip-up clasts. 
The immature pebbles in the top row were collected close to 
their origin; the sub-mature pebbles in the middle row were 
collected further downstream; the smooth, mature pebbles 
in the bottom row were collected on a nearby sandbar. As 
these pebbles eroded, their shapes became rounder, an effect 
quantified in the next figure. The bars indicate, in each row, 
a length of 2 cm. 



in the world after the Bay of Fundy, Canada. During 
the spring tide periods the upper part of the tidal flats 
collects a muddy sediment. This mud dries up during 
the following neap tide period where the sediments are 
exposed to the air. In certain areas, between the large 
cquinoxial neap tides, the exposure of mud sediments 
to air may last for several months. During this period, 
this mud layer will develop a vast network of desiccation 
cracks. This network then leads to fragmented plates of a 
polygonal shape with 20 to 40 cm size. During the next 
spring tide period, the plates in the erosional area can 
be eroded by tide currents, thus re-incorporated into the 
sedimentary cycle. During the tide process these clasts 
are progressively eroded over many months. Thus the 
mud cohesion allows enough observation time for the life 
of a clast to be observed within a distance of order of one 
hectometer, from the original erosional area down to the 
latest stages of abrasion. 
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FIG. 9: Cumulative curvature distribution, f(K), for the 
average shapes of Mont St.-Michel rip-up clasts. Even the 
roundest shapes remained less circular than the final shape in 
the laboratory study. 



We have analyzed the shapes of three classes of rip- 
up clast photographed at three distinct locations on the 
tidal flats near Mont St.-Michel. The first is large sub- 
angular cobble found near the site of formation. Ten 
samples were examined; for these immature pebbles, the 
average perimeter is 550 mm. The second class is medium 
sub-mature pebbles found further "downstream" . As a 
result of erosion, these pebbles are smaller and smoother 
than the cobble. Thirty five samples were examined; for 
these, the average perimeter is 180 mm. The third class 
is rounded, mature pebbles found on a nearby sand bar. 
The relation of these pebbles to the other two classes is 
not clear. Seventeen samples were examined; for these, 
the average perimeter is 220 mm. 

Typical photographs for each of these classes are shown 
in Fig. [8] The average of the cumulative curvature distri- 
bution for all samples in each class is shown in Fig-EI The 
angularity of the large cobble is reflected in the breadth 
of the curvature distribution. Roughly a quarter of the 
perimeter has negative curvature, and roughly a tenth 
has curvature five times greater than the average. For 
the other two classes, the curvature distribution is pro- 
gressively more narrow. The relative steepness of the 
/(-ftT)'s shows that all of these shapes are less round than 
the final pebbles produced by the laboratory erosion ma- 
chine. 

The width of the curvature distribution can be speci- 
fied quantitatively by the standard deviation, a. Results 
are normalized by the average curvature, and are shown 
for the field and laboratory pebbles in Table [I] The peb- 
bles with steeper f(K) indeed have smaller widths. For 
example, the width for the immature field cobble is about 
3-4 times that of the average laboratory pebble. While 
the dimensionless width of the distribution, a/(K), is a 
useful number for comparisons, it does not distinguish 
between curvature distributions of different shape. The 
actual functional form of the curvature distribution can 
be specified to some extent by comparing its moments 
with that of a Gaussian. In particular, the "skewness" 
and "kurtosis" are dimensionless numbers defined by the 
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third and fourth moments, respectively, in such a way as 
to vanish for a perfect Gaussian. The results in Table Q] 
show that the four classes of pebbles have curvature dis- 
tributions of four distinct forms. Of these, the laboratory 
pebbles are closest to a Gaussian. 



Class 


Perimeter (mm) 


a/(K) 


Skewness 


Kurtosis 


immature 


550±100 


2.4±0.3 


1.4±0.3 


3.8±1.4 


sub-mature 


180±60 


2.0±0.6 


0.0±0.9 


3.3±2.9 


mature 


220±70 


1.4±0.3 


-0.2±0.5 


1.3±1.0 


lab- final 


122±25 


0.8±0.1 


0.1±0.5 


1.0±1.5 



TABLE I: Characteristics of curvature distribution for the 
three classes of field pebbles. Final laboratory pebble shape 
is added for comparison. 



V. CONCLUSION 

We have studied the formation of two-dimensional peb- 
ble shapes. As in Ref. we introduced a local descrip- 
tion of the erosion process, based on the distribution 
function of the curvature, measured along the pebble con- 
tours. This description captures both the local character 
of the erosion events, and the statistical nature of the 
erosion process. 

For pebbles generated in the laboratory, we have shown 
that the curvature distribution has two important prop- 
erties. First, the erosion drives the distribution towards 
a stationary form. When this stationary state is reached, 
the pebble contour still changes but, within small fluc- 
tuations, its curvature distribution remains the same, 
provided that the curvature is normalized by its average 
value. Secondly, we have found that the final stationary 
form of the distribution is independent of the original 
pebble shapes. This not only shows that the curvature 
distribution is a property of the erosion process itself, but 
it also opens the interesting possibility of establishing a 
classification of different erosion processes according to 
the type of curvature distribution they generate. 

For pebbles collected in the field, we have made a first 
attempt to study a special class of rip-up clasts from 
the St. -Michel bay. These mud pebbles can be collected 
at very different erosion stages within a relatively small 
area of the tidal flats. We showed that the curvature 
distribution sharpens with the wearing degree, without 
getting however as sharp as the distribution obtained in 
the laboratory experiments. 

The results presented in this experimental paper sug- 
gest a number of directions for modeling the formation of 
flat pebbles. Of central importance is the intrinsic statis- 
tical nature of the erosion process itself. As first hinted 
in Ref. a sequence of cuts of a noiseless, deterministic 
nature typically leads to a trivial curvature distribution 
like that of a circle. We also demonstrated in Ref. 8] that 
a "cutting" simulation, with an appropriate distribution 



of cutting lengths, acting most strongly on regions of high 
curvature in accord with Aristotle's intuition can re- 
produce the curvature distribution from the laboratory 
experiments. We will address these and other questions 
relevant for the theoretical modeling of pebble formation 
in a forthcoming paper. 
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APPENDIX: CURVATURE ANALYSIS 

The goal of this section is to provide a detailed, prac- 
tical description of two means to measure the local cur- 
vature at each point along the pebble contour. In both 
cases the starting point is an image of the pebble. For 
our work we use a digital camera Canon Power Shot Gl 
with a resolution of 1024 x 768 pixels. It should be just 
as effective to scan conventional photographs, or even to 
scan pebbles themselves. To determine the {x,y} coor- 
dinates of the contour, we import the images into NIH 
Image which includes routines for finding edges and 
for skeletonizing the result. An example of the digitized 
pebble contour, and the smooth reconstructions to be dis- 
cussed below, is given in Fig. ^| We show the contour 
in pixel units, since the actual calibration is needed only 
to determine size, not shape. Note from the inset that 
the contour points are indeed pixelized and skeletonized, 
with each point having only two neighbors located at ei- 
ther ±1 or units away in the x- and y-directions. If the 
digitization process is faithful, then the uncertainty in 
each pixel is about ±0.5 units in each direction. This is 
not small compared to the distance between neighboring 
pixels, so a smoothing or fitting routine is necessary to 
reconstruct the actual contour and thereby extract the 
curvature distribution. 

To illustrate the difficulties of extracting curvature, let 
us begin with two methods that are, in fact, unsatis- 
factory. It is perhaps tempting to simply smooth the 
data, replacing each point with a weighted average of 
neighbors lying within some window. Weights could be 
cleverly chosen to de-emphasize points at the edge of the 
window, for example. This fails, however, since it's far 
from obvious how to choose a suitable window size. For 
instance, the pixelized representation of the straight sec- 
tion given by y = O.lx for < x < 10 is a step function 
y = 0(1) for x < (>)5. A large window would be needed 
to even approximately reconstruct the original line. How- 
ever, such a sufficiently large smoothing window would 
erase fine features if applied elsewhere along the contour. 
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FIG. 10: Reconstruction of a smooth pebble contour from 
the pixelized digital representation. The solid curve is based 
on fitting to cubic polynomials at each point. The dashed 
curve is based on the iteration scheme of Fig. 1111 



Since smoothing niters provide no feedback on quality, vi- 
sual inspection of the result would be necessary to choose 
an optimal window size at each point along the contour. 
This is not only subjective, but rather impractical. As 
an alternative, it is perhaps tempting to implement an 
automated version of Wentworth's curvature gauge 
This is a device with circular notches of various diam- 
eters into which portions of a pebble may be pressed. 
The computational analogue would be to find the best 
nonlinear-least-squares fit to a circle at each point along 
the pebble contour. As with smoothing, one difficulty is 
to find the optimal window over which to do the fitting. 
A compounding difficulty is that essentially nowhere is 
the pebble exactly circular, so even with 'the' optimal 
window there is substantial disagreement between data 
and fit. A spectacular example of this problem is at an 
inflection point, where the curvature changes from posi- 
tive to negative. 

To overcome such difficulties, we propose to fit digi- 
tized contour data to a cubic polynomial at each point 
along the contour. A cubic is the lowest order polynomial 
needed in order to avoid systematic error when the cur- 
vature varies gradually across the fitting window, which 
is the usual case. In order to avoid having to rotate 
the coordinate system to ensure that the contour y(x) 
is a single-valued function, we instead convert to polar 
coordinates. Thus we define the origin by the center- 
of-mass of the contour and perform fits to r(6) where 
r = \J x 2 + y 2 and 9 = tan _1 (y/x). Once a satisfactory 
fit is achieved, the curvature may be deduced from the 
value and derivatives of the cubic polynomial by 

K= (r2+rg)3/2 ' (A - 1} 



where rg = dr/d6 and rgg = d 2 r/d9 2 [29| . 

Two tricks seem necessary to achieve satisfactory re- 
sults. The first is to weight the data most heavily near 
the center of the window. We use a Gaussian weighting 
function with a standard deviation equal to 1/4 of the 
width of the window. This ensures that points at the 
edges have essentially no influence. Therefore, the fitting 
results do not vary rapidly as the window is slid along the 
contour. This guarantees that the reconstructed curve 
and its first two derivatives are continuous, which is a 
crucial requirement for measuring the curvature. 

The second trick is to choose the window size appro- 
priately. This is actually the most difficult and subtle 
aspect of the whole problem. If the window is too small, 
then the fit will reproduce the bumps and wiggles of the 
pixelization process; usually the curvature will be over- 
estimated. If the window is too large, then the fit will 
significantly deviate from the data; usually the curvature 
will be underestimated. And while the curvature tends 
to decrease systematically with window size, there is in 
general, unfortunately, no plateau between these two ex- 
tremes where the curvature is relatively independent of 
window size and hence clearly represents the true value. 
To pick the window appropriately requires careful un- 
derstanding of the numerical fitting procedure and the 
feedback it provides. Since the fitting function is a poly- 
nomial, the minimization of the \ 2 total square deviation 
from the data reduces to solving a set of linear equations. 
This in turn reduces to inverting a matrix. If the window 
is too small then the fit will be 'ambiguous' in the sense 
that x 2 is small but the error in the fitting parameters is 
large. Mathematically, the matrix to be inverted is essen- 
tially singular. A good strategy is therefore to start with 
a small window and increase its size until the matrix is 
no longer singular. This can be accomplished using a lin- 
ear least-squares fitting routine based on singular-value 
decomposition 30]. However, the uncertainty in fitting 
parameters for the first suitable window is generally too 
large (nearly 100%). So we increase the window size two 
pixels at a time until the error in curvature has been re- 
duced by a factor of ten. This defines the largest suitable 
window, beyond which systematic errors due to incorrect 
functional form begin to appear. We have not been able 
to define the largest suitable window based on the value 
of x 2 . For the final result, we take a weighted average of 
the cubic fit parameters over all suitable windows, where 
the weights are set by the uncertainties in fitting param- 
eters as returned by the fitting routine. An example of 
a reconstructed contour from this procedure is shown by 
the solid curve in Fig. ^| The recontruction is satis- 
factorily smooth; also, it clearly avoids the pixel noise 
without smoothing over significant small-scale features 
in the contour. 

Since the cubic fitting method is rather involved, and 
since the choice of window sizes is still slightly subjective, 
we have developed an alternative method. The starting 
point is the fact that the actual digital representation 
of the contour depends on the location and orientation 
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FIG. 11: Trial integrated curvature distributions from the 
iterative smoothing scheme, vs the number of points kept. 
Keeping either 4 or 5 points seems optimal: the resulting 
distributions are identical and they agree with that based on 
cubic fits. 

of the pebble with respect to the grid of pixels. If the 
pebble were shifted or rotated, then the pixelized repre- 
sentation would be slightly different. For example, imag- 
ine the pixelization of a line making various angles with 
the grid. Perhaps the ideal experimental measurement 
procedure would be to systematically reposition the peb- 
ble, pixelize, then compute the average of all such rep- 
resentations. However this procedure does not lend it- 
self to automatization, and would be impractically time- 
consuming. Instead, we propose to do more or less the 
same thing numerically. The idea is to take the current 



best guess for the contour, pixelize it with respect to a 
random grid position and orientation, then use the new 
representation to update the best guess. When iterated, 
this procedure converges to a satisfactory reconstruction 
of the actual contour with two provisos. First, at each 
step, we locally smooth the trial pixelization by replacing 
each point by its average with its two immediate neigh- 
bors. Second, we keep only every fourth or fifth point 
in the original pixelized data and perform all operations 
on this subset. When done, we compute the curvature 
literally by the change of slope with respect to arclength 
using the straight segments between adjacent points. 

The cumulative curvature distribution given by this it- 
eration scheme is shown in Fig. ^] as a function of the 
number of points kept. When too many are kept, the re- 
constructed curve follows the bumps and wiggles of the 
original pixelization too closely; the curvature distribu- 
tion is too broad. When too few are kept, the recon- 
structed curve incorrectly smooths over small-scale fea- 
tures; the curvature distribution is too narrow. When 
only every fourth or fifth point are kept, the distributions 
are nearly equal; furthermore, they are indistinguishable 
from that given by the cubic polynomial fitting. The ac- 
tual reconstructed contour is also shown in Fig. 1101 The 
plateau in the curvature distribution vs number of points 
kept, and the good agreement with the other method, 
both give confidence in this new iterative reconstruction 
scheme. 
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